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In this work we compare the results of the Gross-Pitaevskii and modified Gross-Pitaevskii equa- 
tions with ab initio variational Monte Carlo calculations for Bose-Einstein condensates of atoms 
in axially symmetric traps. We examine both the ground state and excited states having a vortex 
line along the z-axis at high values of the gas parameter and demonstrate an excellent agreement 
between the modified Gross-Pitaevskii and ab initio Monte Carlo methods, both for the ground and 
vortex states. 
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I. INTRODUCTION 

Most theoretical studies of Bose-Einstein condensates (BEC) in gases of alkali atoms confined in magnetic or optical 
traps have been conducted in the framework of the Gross-Pitaevskii (GP) equation 0. The key point for the validity 
of this description is the dilute condition of these systems, i.e., the average distance between the atoms is much larger 
than the range of the inter-atomic interaction. In this situation the physics is dominated by two-body collisions, well 
described in terms of the s-wave scattering length a. The crucial parameter defining the condition for diluteness is 
the gas parameter x(r) = n(r)a 3 , where n(r) is the local density of the system. For low values of the average gas 
parameter x av < 10~ 3 , the mean field Gross-Pitaevskii equation does an excellent job (see for example Ref. [2| for a 
review). However, in recent experiments, the local gas parameter may well exceed this value due to the possibility of 
tuning the scattering length in the presence of a Feshbach resonance 0, 0] • 

Under such circumstances it is unavoidable to test the accuracy of the GP equation by performing microscopic 
calculations. If we consider cases where the gas parameter has been driven to a region were one can still have a 
universal regime, i.e., that the specific shape of the potential is unimportant, we may attempt to describe the system 
as dilute hard spheres whose diameter coincides with the scattering length. However, the value of x is such that the 
calculation of the energy of the uniform hard-sphere Bose gas would require to take into account the second term in 
the low-density expansion |5j of the energy density 



E _ 2Trn 2 ah 2 
V m 



128 



na 

7T 



1/2 



(1) 



where m is the mass of the atoms treated as hard spheres. For the case of uniform systems, the validity of this 
expansion has been carefully studied using Diffusion Monte Carlo and Hyper-Netted-Chain techniques Q- 

The energy functional associated with the GP theory is obtained within the framework of the local-density approx- 
imation (LDA) by keeping only the first term in the low-density expansion of Eq. Q 
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where 



Vtrap(r) = 2 m ( u x x + U) ±.V + UJ z z ) 
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is the confining potential defined by the two angular frequencies u>± and u z . The condensate wave function W is 
normalized to the total number of particles. 

By performing a functional variation of Eqp [if] with respect to \P* one finds the corresponding Euler-Lagrange 
equation, known as the Gross-Pitaevskii (GP) equation 



— V 2 + K rap (r) + 
2m 



m 



(4) 



2 



where ijl is the chemical potential, which accounts for the conservation of the number of particles. Within the LDA 
framework, the next step is to include into the energy functional of Eq. J5J) the next term of the low density expansion 
of Eq. (JJJ. The functional variation gives then rise to the so-called modified GP equation (MGP) Q 
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The MGP corrections have been estimated in Ref. || in a cylindrical condensate in the range of the scattering lengths 
and trap parameters from the first JILA experiments with Feshbach resonances. These experiments took advantage 
of the presence of a Feshbach resonance in the collision of two 85 Rb atoms to tune their scattering length 

Fully 

microscopic calculations using a hard-spheres interaction have also been performed in the framework of Variational 
and Diffusion Monte Carlo methods El El El 

In this work we compare the results of the GP and MGP equations discussed above, Eqs. (@J and J5J, with variational 
Monte Carlo (VMC) calculations for axially symmetric traps in a region [x > 1CP 3 ) where the validity of the GP 
equation is not clear. We examine both the ground state and excited states having a vortex line along the z-axis. 

In the next Section we present our numerical approaches together with a discussion of ground state properties. In 
Sect. Hill we proceed to study several trial wave functions to describe the excited state with one vortex. A comparison 
between VMC and the GP and MGP equations is done. We summarize our results in Sect. II VI 



II. NUMERICAL APPROACHES AND GROUND STATE PROPERTIES 



The starting point of the Monte Carlo calculations is the Hamiltonian for N trapped interacting atoms given by 
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The two-body interaction Vi n t(| r$ — Yj |) between the atoms is described by a hard-core potential of radius a, where a 
is the scattering length. The atoms are thus treated as hard spheres. The next step is to define a trial wave function, 
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where F(l, N) is a many-body correlation operator applied to the mean-field wave function $mf- The advantage 
of using a correlated trial wave function lies in the fact that non-perturbative effects, as the short-range repulsion 
between atoms may be directly incorporated into the trial wave function. The simplest correlation operator has the 
Jastrow form [13, 
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In our variational calculations we use a two-body correlation function which is the solution of the Schrodinger equation 
for a pair of atoms at very low energy interacting via a hard-core potential of diameter a. The ansatz for the correlation 
function f(r) reads 



f(r) = 
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This type of correlation, besides being physically motivated, has been successfully used in Refs. [Tot ITll| to study both 
spherically symmetric and deformed traps. These authors have also explored the quality of this correlation function by 
comparing Variational Monte Carlo (VMC) and Diffusion Monte Carlo (DMC) calculations for the case of spherically 
symmetric traps E3|) with a good agreement between the VMC and DMC results. 

The deformation of the trap is incorporated in the mean field wave function ^mf, which is taken as the product 
of N single particle wave functions 
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where a is taken as the variational parameter of the calculation, and A{a) — (a/n) 3 / 4 is the normalization constant. 
The parameter A = lo z /uj± is kept fixed and set equal to the asymmetry of the trap. In this way the mean- field wave 
function $mf has all the particles in the condensate, the latter being described by the wave function 99. 
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The evaluation of the expectation value of the Hamiltonian with this correlated trial wave function provides an 
upper bound to the ground state energy of the system 

(jgr | H | gg) 
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This expectation value has been evaluated by the Metropolis Monte Carlo method of integration [TEl Il6| . 

The energy obtained with the Hamiltonian of Eq. @ can be directly compared with the output of the GP and 
MGP equations, see Eqs. J3J and JSJ. The Gross-Pitaevskii equations represent a mean-field description, with all the 
atoms in the condensate. In fact, the additional correlations which are taken into account in the second order term 
of the low-density expansion of the energy, see Eq. 1|T|1. are incorporated in the density functional and, therefore, in 
the solution of the MGP equation. In contrast, the Monte Carlo calculation incorporates explicitly the inter-atomic 
correlations and therefore one could in principle find the natural orbits and extract the occupation of the condensate 

M- _ 

The GP and the MGP equations have been solved by the steepest descent method |l7J for the deformed harmonic 
oscillator trap previously described in Eq. @. An initial deformed trial state is projected onto the minimum of the 
functional by propagating it in imaginary time. In practice, one chooses a small time step At and iterates the equation 

^(r,t + At)^^(r,t)-AtH^(r,t) (12) 

by normalizing \I> at each iteration. When the gas parameter becomes large, the time step which governs the rate 
of convergence should be taken accordingly small. Convergence is reached when the chemical potential becomes a 
constant independent of the position, see Eqs. (QJ and (JSJ. 

For the comparison of the results obtained with the different GP type equations and the variational Monte Carlo 
calculations, we consider a disk-shaped trap with A = uj z /uj± = see Ref. 18j. We have fixed the scattering length 
to a = 35aRb, with aRb = lOOao, ao being the Bohr radius. We set the number of confined atoms to N = 500 in order 
to keep the amount of computing time acceptable when using the Monte Carlo method. All the numerical results are 
given in units of the harmonic oscillator length a± = (h/mcj^) 1 ^ 2 and the harmonic oscillator energy huj±. 

First we analyze the GP and the MGP results reported in Table|U For a scattering length a = 35etRb, the corrections 
of the MGP approach to the chemical potential are of the order of 20%. The energy corrections are also relevant and 
it is interesting to study the different contributions to the energy. The kinetic energy is given by 

Skin = ^ / dr | V*(r) | 2 , (13) 
while the harmonic oscillator energy due to the trapping potential reads 

Eko = jJ dr (,1(, 2 + y 2 ) + u 2 z z 2 ) | *(r) | 2 , (14) 
and the interaction energies E\ and E2 are given by 
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The virial theorem is used to establish a relation between the different contributions to the energy, viz., 

2E kin - 2E HO + 3£a + \e 2 = 0, (17) 

which serves as a proof of the numerical accuracy of the solution of the GP equations. The results in Table [I] show 
that this test is well satisfied by all calculations. 

Notice that the kinetic energy associated with the mean field descriptions is not negligible, indicating that the 
regime where the Thomas-Fermi approximation to the GP equation is valid has not been reached. In this limit, the 
chemical potential is 



= i(15aA^A) 2/5 ^ (18) 
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where a = a/a± is the dimensionless scattering length, and the energy per particle Etf/N — 5/itf/7. In this 
approach we have E^f/N — 9.03Tilu± and /j,tf = 12.64Juj±. Both these values differ from the values reported in 
Table [I] However, this approximation can still be used to estimate the peak value of the gas parameter, namely 

4 fe = n (0)a 3 = -^(15aiVA) 2 / 5 a 2 , (19) 

which yields x^ F = 0.023. At this rather large value of the diluteness parameter, the corrections brought by the MGP 
equation to the GP results are expected to be relevant However, x is low enough to allow for a mean- field 

approach (as it is the case of the MGP equation). For such density regimes, a mean- field approach provides a rather 
good description when compared with a microscopic calculation j3| ■ 

The variational Monte Carlo results are also given in Table [I] and they show a close agreement with the results 
provided by the MGP equation. Notice that in this approach, and using the Hamiltonian of Eq. (JBJ, the potential 
energy is zero since the wave function is strictly zero inside the core. The total energy in this case is distributed 
between -Eho an d the true kinetic energy. Actually the only energies that can be directly compared with the GP 
results are the total and the harmonic oscillator energies. 

The Monte Carlo results obtained with the Metropolis algorithm take into account the energy of 27000 configura- 
tions, grouped in 90 blocks of 300 movements. At each Monte Carlo step we move all the particles and the acceptance 
is around 58 %. A thermalization process is incorporated at the beginning of the Monte Carlo process and before each 
block. In the Monte Carlo calculation we have used the Pandharipande-Bethe prescription for the kinetic energy [T^| . 
which produces a smaller variance. To a give a feeling for the numerical accuracy of our VMC results, we list here 
GP, MGP and VMC results in the dilute limit. We employ N = 500 particles and a scattering length for 87 Rb 
considered by Dalfovo and Stringari in Ref. [Tflj , which in units of the oscillator parameter perpendicular to the 2-axis 
is 4.33 x 10~ 3 cu. We obtain energies in units of the oscillator energy of 3.3032, 3.3080 and 3.3242(1) for GP, MGP 
and VMC calculations, respectively. The VMC results are for an optimum variational parameter a = 0.475. Taking 
into account that the two-body correlation has been kept fixed, and that the only variational parameter is a, these 
results indicates that our ansatz for the variational wave function is a viable one. Actually, as the reader will notice 
from the discussion below, this discrepancy of roughly 0.5% is of the same relative order as for the higher density 
cases reported here. 

In the minimization process we keep fixed the parameter A in the single particle wave function of Eq. I |10|l . i.e., 
we assume that the deformation of the trap is transferred to the wave function, and vary only a. At the minimum, 
a = 0.7687. One can also explore the effects of the correlations in the density profiles. These profiles, which represent 
a column density defined according to 

n c (r±) = Jdzn(r±,z) (20) 

and normalized such that 2tt J dr± r± n c (r±) = 1, are shown in Fig. ^ for the various approximations used in this 
work. The repulsive character of the correction term of the MGP equation translates into a decrease of the value 
of the column density at the origin and an increase of the size of the condensate 0, || . This gives a slightly more 
extended profile for the MGP approach compared with both the GP and the VMC results. As one can see from Fig.^ 
there is a much better agreement between the Monte Carlo and MGP profiles than with the corresponding profile 
from the GP calculation, particularly at small values of the radial distance where the density is larger. 

The good agreement between VMC and MGP does not guarantee that these methods give a good description of 
the system. However, as it was shown in Ref. [Tl| for the case of spherical traps, the improvements introduced in 
the trial wave function by a Diffusion Monte Carlo calculation, which in principle allows for an exact solution of the 
many-body problem, are rather small and the variational wave function of Eq. i|10|) provides a very good description 
of the system. Therefore we assume that the same will be true for deformed traps. Furthermore, for these values 
of the diluteness parameter, the MGP equation is very useful to calculate the energy, chemical potential and the 
density profiles of the ground state of the system for condensates with larger number of particles, which would be 
computationally prohibitive for a Monte Carlo calculation. 



III. VORTEX STATES 



The existence of these excited condensate states is crucial to studies of the superfluid behavior of trapped atomic 
condensates. In this section we study the effects of correlations in vortex states. We consider a singly quantized 
vortex line along the z-axis. This means that all the atoms rotate around z-axis with a quantized circulation. The 
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GP equation can easily be generalized to describe this kind of vortex states [2j by using the following ansatz for the 
condensate wave function 

*(r) = ip(r) axp (21) 
where <j> is the angle around the z-axis and k is an integer. This vortex state has a tangential velocity 

u<# = K (22) 

mr± 

where rj_ = \J x 2 + y 2 is the distance to the symmetry axis of the vortex. The number k represents the quantum 
of circulation, and the total angular momentum along the z axis is given by Nk%. Introducing the wave function of 
Eq. (|21|l in the GP energy functional of Eq. © one gets the corresponding GP energy functional for the vortex state 



-BGP+vor[*] = J dr 
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which incorporates a centrifugal term in the density functional, arising from the quantized flow of atoms around the 
vortex core. This term defines a rotational energy 

£rot = £/dr^(r)| 2 . (24) 



The corresponding nonlinear Schrodingcr equation obtained by functional variation is 

ip = [iip. (25) 
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Adding E% to the density functional and after performing a functional variation one gets the corresponding MGP 
equation for the vortex state. 

Based on the virial theorem, one can again derive a relation between the different contributions to the energy: 

2E kin - 2E HO + m + ~£ 2 + 2E Iot = 0. (26) 

The thermodynamic critical angular frequency Q c required to produce a vortex of vorticity k is obtained by com- 
paring the energy of the system in the rotating frame with and without the vortex |20| : 

n c =^[K-£ 0] . (27) 

A main feature of a vortex state is the hole (core of the vortex) that appears in the center of the density profile 
along the rotation axis. From Eq. I|25|) . it is clear that the solution of this equation has to vanish on the z-axis due 
to the presence of the centrifugal term. The size of the core is characterized by the healing length. 

For the microscopic description of the vortex state we use an Onsager-Feynman type trial wave function pH 

<M1, N) = e iK ^ ^ J] /(r^-)*o(l, N), (28) 

3 

where N) is the ground state wave function. The phase factor K^j^i depends on the angular variables 

of the particles and is the equivalent to the phase factor introduced in the mean field description of Eq. 121fl . The 
function f(r±) modulates the density as a function of the radial coordinate r±. We examine three types of f(r±). In 
the first ansatz we use the simple option 

/i(rx)=r±. (29) 

In the second case we consider, 



h{r±) = l-exp(-rj_/d), 



(30) 
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where d is a variational parameter. Notice that for d = 1, the behavior of /2(t\l) for small r± coincides with the 
behavior of f%(r±). Finally the third function, is that of Ref. [23 which has been used in the context of quantum 
liquids, 

/ 3 (rx) = 1 - exp (-(rx/d) 2 ) , (31) 

where d is again a variational parameter. 

These three trial wave functions describe a singly quantized vortex state (k = 1), whose axis lies in the z-direction 
and with a tangential velocity field — h/mr±. The evaluation of the expectation value of the Hamiltonian (Eq. ©) 
with these wave functions is equivalent to calculate the mean value of the Hamiltonian 

h 2 N N ^2 2 N N 

ff = ^E^ 2 +E^+E w*o + E mi * - *i i). (32) 

i—1 i—1 -Lj^ i—1 i<j 

with 9(1, ...N) = Y[j f( r ±.j)9o(l, N). In this way the rotational contribution to the energy has been directly 
incorporated in the Hamiltonian. Minimizing this new problem provides the best energy and wave functions inside 
this subspace of wave functions. In the context of liquid 4 He there have been attempts to perform a full minimization 
allowing for a more general phase function. The analysis indicates that the present procedure provides very accurate 
results |23|. 

We start by discussing the GP and MGP results (obtained by the steepest descent method as done for the 
ground state as well) with an initial condensate wave function 

1>(t) ex /i(r±)* (r). (33) 

It is worth mentioning, as a check of the numerical procedure, that starting with /a(r_L) or fs(r±) to modulate the 
condensate wave function we converge to the same results as with fi(r±). 

As expected, the presence of the vortex increases the chemical potential. Also £ho has a small increase, related to 
the enlargement of the profile due to the presence of the vortex hole. These results are listed in Table[H] Although the 
MGP corrections to the energy are sizable and of the same order as those in the ground state, the critical frequency, 
Ogp = 0.29uj± , is barely affected as both energies, the energy of the vortex state and the ground state energy, are 
shifted by similar amounts, yielding VLmgp = 0.24u^. 

The GP and MGP profiles for the vortex state are shown in Fig. [2] As a consequence of the repulsive character 
of the MGP corrections, the central density of the GP ground-state density profile is higher than the MGP one and, 
therefore, the depth of the hole around the z axis is larger in the GP approach. However, the healing length is almost 
the same. 

As can be seen from Table ITTT1 the Monte Carlo results for the energies are in good agreement with the MGP ones 
for all the trial wave functions considered. This table shows two types of calculations. In the first three rows we list 
the energies obtained by keeping 9q equal to the ground state wave function and performing the minimization with 
respect to the parameter d in the modulating function, except in the case of /i, which has no variational parameters. 
In the second set of results, we perform a minimization allowing to vary also the harmonic oscillator parameter a of 
the wave function \&o- The changes in a and d do not yield significant changes in the computed energy. 

The density profiles seem to be more sensitive to the modulating function, as one can see from Fig. |2 . These 
profiles correspond to the case where the ground state wave function 9q is kept fixed when we minimize the energy of 
the vortex state. For we obtain a radial structure which is not present in the mean field approach |22 • The MGP 
profile shows a broader surface region than the VMC profiles. In the core of the vortex, the MGP profile looks very 
similar to the VMC results with the modulating function /a of Eq. J3UJ). These two results exhibit a smaller healing 
length than the VMC calculation which employs /i . 

From the variational point of view, the best description of the vortex should correspond to the wave function that 
provides the minimum energy. According to this criterion, this corresponds to the trial wave function built with the 
modulating function fi of Eq. I|29|) . 

Finally, in Fig. [31 we plot the density profiles for all VMC calculations, with and without vortices. We note that 
they all provide a similar healing length and that the asymptotic behavior is almost equal for both the ground state 
and the vortex states. 



IV. CONCLUSIONS 

We have compared the results of the Gross-Pitaevskii (GP) and the modified Gross-Pitaevskii (MGP) equations 
with ab initio variational Monte Carlo calculations for Bose-Einstein condensates of atoms in deformed traps. We 
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have studied both the ground-state and excited states having a vortex line along the z-axis. The interatomic potential 
has been characterized by a hard-sphere potential with a radius which coincides with the scattering length used in 
the GP and MGP equations. 

We have performed the calculations for 500 particles. The parameters characterizing the trap and the scattering 
length have been chosen to reach values of the gas parameter where the MGP calculations provide corrections of 
the order of 20% compared with the GP results. It is indeed very interesting that even at such values of the gas 
parameter one can still describe the system in terms of mean-field approaches. We find for example an excellent 
agreement between the MGP and the VMC results, especially for the energies of the ground state and the vortex 
states. The MGP and VMC density profiles for the ground state are also in good agreement. The situation is different 
for the vortex state. Three different trial wave functions produce similar energies but slightly different profiles. In the 
core of the vortex, the MGP profile is close to the profiles obtained with the ansatzes /i and of Eqs. i|29[l and (I3U|> . 
respectively. The latter two wave functions yield also the lowest energies. Whether a Diffusion Monte Carlo (DMC) 
calculation will show a similar trend remains to see. We are planning DMC studies of the systems discussed here. 
Our preliminary DMC calculations for the energy of the ground state show little change with respect to the VMC 
results and, hence, a very good agreement with the MGP results. 

In summary, we would like to point out that the good agreement between the VMC and MGP is rather encouraging, 
and allows for further MGP explorations of vortex states in condensates with both a larger number of interacting 
atoms and large scattering lengths. 
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FIG. 1: Ground state column density n c (rj_) as a function of the distance to the z axis, for N = 500 particles, comparing the 
GP (solid line) and MGP (dashed line) results for a = 35dRb = 0.15155ax. Also shown are the results of variational Monte 
Carlo calculations (line with symbols) . The deformation A = VS and the oscillator lengths are defined as in Refs. 0, 0| . The 
radial distance is given in units of aj_ = (h/muj^) 1 ^ 2 . The column density is dimensionless. See text for further details. 
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FIG. 2: Vortex column density n c (r±) as a function of the distance to the z axis, for N = 500 particles, comparing GP 
(solid line) and MGP (dashed line) results for a — 35aRb = 0.15155ax. Also shown are the results of variational Monte Carlo 
calculations (lines with symbols) using the different Onsager-Feynman ansatzes. The deformation A = v» and the oscillator 
lengths are defined as in Refs. [Tsl Il9( . The radial distance is given in units of a± = (li/mwi) 1 ' 2 . The column density is 
dimensionless. See text for further details. 
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FIG. 3: Column density n c (r±) as a function of the distance to the z axis, comparing the VMC profiles for the vortex state 
corresponding to the different Onsager-Feynman ansatzes (lines with symbols as in Fig.(|2J and the ground state (dashed line 
with full circles). The trap parameters and the scattering length are the same as in the two preceeding figures. The radial 
distance is given in units of a±_ = (h/mu;±) 1 ^ 2 . The column density is dimensionless. See text for further details. 
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TABLE I: Chemical potential and energies in units of hui± from the GP, MGP and VMC calculations for the ground state. 
The scattering length is a = 35a Rb = 0.15155a ± , A = \/8, N — 500. 

M E/N E kin /N Euo/N Ei/N E 2 /N 

GP 12.980 9.496836 0.39495 5.61911 3.4827765 

MGP 15.453 11.06108 0.35353 6.94092 2.516691 1.249938 

VMC 11.12109(14) 4.21520(24) 6.90590(19) 



TABLE II: Chemical potential and energies in units of hu± from the GP and MGP calculations for the one-vortex state with 
the vortex line along the z-axis. The scattering length is a = 35a Rb = 0.15155a_L, A = \/8, N = 500. 







E/N 


E Un /N 


Eho/N 


Ei/N 


E 2 /N 


E rot /N 


GP-lv 


13.187 


9.7835936 


0.42508 


5.74271 


3.403871 




0.21193 


MGP-lv 


15.623 


11.305 


0.37692 


7.03774 


2.482418 


1.223280 


0.18492 



TABLE III: Variational Monte Carlo results obtained with different Onsager-Feynman ansatzes. The results labeled VMC[fl], 
VMC[f2] and VMC[f3] stand for the modulating wave functions in Eqs. I28H . 130H . and I31H respectively. 





Q 


d 


E/N 


E kin /N 


Eho/N 


E ro t 


VMC[fl] 


0.7685 




11.33432(18) 


4.14804(26) 


7.02175(19) 


0.164527(44) 


VMC[f2] 


0.7685 


1.175 


11.36273(18) 


4.23594(31) 


6.93987(24) 


0.186912(58) 


VMC[f3] 


0.7685 


0.425 


11.39171(18) 


4.28845(30) 


6.91368(23) 


0.189580(30) 


VMC[fl] 


0.775 




11.33415(17) 


4.18634(29) 


6.98213(22) 


0.165679(60) 


VMC[£2] 


0.745 


1.425 


11.35457(15) 


4.07816(31) 


7.09696(25) 


0.179446(93) 


VMC[f3] 


0.745 


0.550 


11.38683(19) 


4.14902(33) 


7.06446(26) 


0.173350(26) 



